First principles theory of the EPR g-tensor in solids: defects in quartz 
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A theory for the reliable prediction of the EPR g-tensor for paramagnetic defects in solids is pre- 
sented. It is based on density functional theory and on the gauge including projector augmented 
wave (GIPAW) approach to the calculation of all-electron magnetic response. The method is val- 
idated by comparison with existing quantum chemical and experimental data for a selection of 
diatomic radicals. We then perform the first prediction of EPR g-tensors in the solid state and 
find the results to be in excellent agreement with experiment for the E[ and substitutional P defect 
centers in quartz. 
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Electron paramagnetic resonance (EPR), also known 
as electron spin resonance (ESR), is the most powerful 
spectroscopic technique for the study of paramagnetic 
defects in solids. Indeed, defect centers are often named 
directly after their EPR spectra. Applications of EPR 
extend to any situation where there are unpaired elec- 
trons, including the understanding of reactions involving 
free radicals in both biological and chemical contexts or 
the study of the structure and spin state of transition 
metal complexes. 

EPR spectra of spin 1/2 centers are made up of two 
contributions: (i) the hyperfine parameters, which can 
be computed from the ground state spin density, and 
have been used to connect theoretical studies of defects 
to available experimental data Jl]"^], and (ii) the g-tensor. 
Only recently have there been attempts to calculate the 
g-tensor in molecules from first principles using den- 
sity functional theory (DFT) [[?]||]. However, these ap- 
proaches are valid only for finite systems and, thus, are 
not useful for the calculation of the g-tensor for param- 
agnetic defects in solids, except possibly within a cluster 
approximation. In the absence of a predictive scheme, 
experimentally determined g-tensors are, of necessity, in- 
terpreted in terms of their symmetry alone, leaving any 
remaining information unexploited. A reliable, first prin- 
ciples approach to the prediction of g-tensors in solids, in 
combination with structural and energetic calculations, 
would access this information, and could be used for 
an unequivocal discrimination between competing micro- 
scopic models proposed for defect centers. In this let- 
ter we describe an approach for the calculation of the 
g-tensor in extended systems, using periodic boundary 
conditions and super-cells. 

In a previous paper || we have shown how to com- 
pute the all-electron magnetic linear response, in finite 
and extended systems, using DFT and pseudopotentials. 



To achieve this we introduced the gauge including pro- 
jector augmented wave (GIPAW) method, which is an 
extension of Blochl's projector augmented wave (PAW) 
method |0J. In Ref. [§, we used GIPAW to compute 
the NMR chemical shifts in molecules and solids. Here, 
we apply the GIPAW approach to the first principles 
prediction of EPR g-tensors for paramagnetic defects in 
solids. We validate our theory and implementation for 
diatomic radicals, for which both all-electron quantum 
chemical calculations and experimental data exist. As, 
until now, there have been no first principles calculations 
of g-tensors in solids, we validate our method in the solid 
state by a direct comparison with experiment. In particu- 
lar, we interpret, from first principles, the EPR spectrum 
of the well characterized and technologically important 
defects in quartz, the E[ and P substitutional centers. 

The g-tensor is an experimentally defined quantity, 
arising from the recognition that the EPR spectrum can 
be modeled using the following effective Hamiltonian, bi- 
linear in the total electron spin S, and the applied uni- 
form magnetic field or nuclear spins, B and I/, respec- 
tively: 



;S ■ g • B + S ' A I ■ h 



(1) 



Here, and in the following, atomic units are used, a is the 
fine structure constant, and the summation I runs over 
the nuclei. The tensors A/ are the hyperfine parame- 
ters (a PAW based theory for its calculation has been 
described elsewhere by Van de Walle and Blochl 
and the tensor g is the EPR g-tensor. 

In order to calculate the g-tensor we start from the 
electronic Hamiltonian which includes terms up to order 
a 3 , in the presence of a constant external magnetic field 
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The summations over i and j run over the electrons and 
Hz, -ffz-KE, -ffso: and -ffsoo are the electron Zeeman, 
the electron Zeeman kinetic energy correction, the spin- 
orbit, and the spin-other-orbit terms respectively: 
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The constant g' is related to g e , the electronic Zeeman 
g-factor in vacuum, by g' = 2(g e — 1), and A(r) = ^Bxr 
is the vector potential. 

Starting from the Hamiltonian of Eq. (0) , we can ex- 
pand the total energy in powers of a, up to 0(a 3 ), using 
perturbation theory. In the resulting expansion, the term 
bilinear in Si and B is identified as the first term of Eq. 
([j]). This term can be rewritten within the formalism of 
spin polarized DFT to obtain an explicit expression for 
the g-tensor: 

g = g c + Ag z _KE + Ag so + Ag S oo = g c + Ag, (4) 
where g c = g c I, I being the identity matrix, and: 



Ag z _ KE = -« 2 g e (Tf -lf)I 



Ag so • B = | g'| d 3 r\jf\r) x V^»(r) 
-jf'WxV^fr)] 



(5) 



(6) 



Agsoo • B = 



2 J d 3 rBW(r)[ P ; o) (r)-pf(r)] (7) 



Here f denotes the majority spin channel and p ® (r) 



T\ \ and V^(r) are the unperturbed electron 
probability-density, kinetic energy, and Kohn-Sham po- 
tential of the j-spin channel, respectively, ji^fr) is the 
electronic charge-current linearly induced by the constant 
magnetic field B in the j-spin channel. Finally, B^^r) is 
the magnetic field produced by the total induced current, 
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[j| (r) +j| 1-l (r)], which we correct for self-interaction by 
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removing the contribution from the current of the un- 
paired electron, [jj^(r) — jj^(r)]. 

We can interpret the physical origin of deviation of the 
g-tensor from its value in vacuum. The spin-other-orbit 
correction, Agsoo, describes the screening of the exter- 
nal field B by the induced electronic currents, as expe- 
rienced by the unpaired electron. The unpaired electron 
itself is not at rest and in the reference frame of the un- 
paired electron the electric field due to the ions and to 
the other electrons is Lorentz transformed so as to ap- 
pear as a magnetic field. The interaction between the 
spin of the unpaired electron and this magnetic field re- 
sults in the the spin-orbit correction, Agso |@- Finally, 
the electron Zeeman kinetic energy correction, Ag z _KE, 
is a purely kinematic relativistic correction. 

Eqs. (|||^) show that the evaluation of the g-tensor 
requires, besides ground state quantities, the linear mag- 
netic response currents j^'(r). Mauri, Pfrommer and 
Louie Q showed how to calculate the magnetic response 
of a system of electrons in an infinite insulating crystal, 
and our recent paper || reformulated this so as to be 
strictly valid for non-local pseudopotentials, and to re- 
produce the valence all-electron currents even within the 
pseudisation core region. An accurate description of the 
all-electron currents in the core regions is essential for the 
evaluation of the SO term, Eq. (0). Indeed, the domi- 
nant contribution to the integral in Eq. (^|) comes from 
the core region as a result of the divergence of (r) at 
the nuclei. 

Using our GIPAW approach to the calculation of all- 
electron magnetic response using pseudopotentials, de- 
scribed in detail in Ref. 0], we break the SO term into 
three parts which derive from the three GIPAW contri- 
butions to the induced current, Eq. (34) of Ref. M: 



Agso = Ag, 
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(8) 



The Agso° term is evaluated from Eq. (^) using a spin 
dependent version of the j^ e (r) of Ref. § and a lo- 
cal Kohn-Sham potential consisting of a sum of the self- 
consistent contribution to the local potential and the lo- 
cal parts of the pseudopotentials. 

The diamagnetic correction term AggQ can be evalu- 
ated from the ground-state valence pseudo wavefunctions 
^ , ^°|) using the following expression : 

Agf d -B= ]T (*%\pi, n K m <Pi, m \&%) 

I,o,n,m 

(*SlPi,n>eU<pj, m |*2>. (9) 

The summation o is over occupied states. The projec- 
tor functions \pr t n) are defined in Ref. || and satisfy 

(jpi,n\4>I',m) = di,i'Sn,m, where \<t>i, n ) are a set of pseudo- 
partial-waves corresponding to the all-electron partial 
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waves \<j>i,n). The projector weights e T nm are given by 
the following atom centered integrals: 



-[(^j,„|(B x r) x WV{r)\4> hm ) 
-(0j, n |(Bxr) x W(r)|&, ro >] (10) 
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The potentials V(r) and V(r) in Eqs. (|L0|) and (|llj) are 
the screened atomic all-electron and local channel pseu- 
dopotentials respectively. 

The evaluation of the paramagnetic correction term 
Ag S Q is more involved as it requires the first order lin- 
ear response wavefunctions. However, the required eval- 
uation can be described by analogy with the calculation 
of paramagnetic correction to the NMR chemical shifts, 
aVjiPAWi replacing the weights f 7 m in Eq. (60) of Ref. 
If by 



f 7 = 



VI, n\ 



lg "L dV(r) 
2 r dr 1 



- g'LdV(r)~ 

h lM^-—^ r \<t>I,m), 
(11) 



where L is the angular momentum operator. 

The electron Zeeman kinetic energy correction term 
Agz-KE is evaluated by combining a straightforward 
PAW correction with the quantity evaluated from the 
ground-state pseudo valence wavefunctions using Eq. (||). 
In this work the SOO term is evaluated from the induced 
field BW(r) derived fr om the bare induced current, and 
the spin density due to the pseudo wavefunctions. It is 
expected that a full GIPAW treatment would result in 
only minor corrections since (i) the SOO term is small 
in comparison to the SO term, and (ii) both the induced 
field and the spin density do not diverge at the nuclei. 

To validate our new expressions for the evaluation of 
the g-tensor, and our implementation of them into a 
parallelized plane-wave pseudopotential code, we com- 
pare with the all-electron gauge including atomic orbital 
(GIAO) DFT results obtained by Schreckenbach and 
Ziegler for a series of diatomic radicals. We use their 
calculated bond lengths for the dimers, but approximate 
the isolated dimers by using large super-cells. Troullier- 
Martins pseudopotentials jl4| and the (spin polarized) 
generalized gradient approximation due to Perdew et al 
|jl5f (GGA-PBE) are used throughout our calculations. 
Table | shows the excellent agreement between our two 
approaches. The exception is the AlO radical, for which 
we obtain much closer agreement with experiment. The 
otherwise close agreement between these two very differ- 
ent approaches suggests a technical rather than funda- 
mental problem in the GIAO calculation for AlO. Com- 
parison with experiment is made through Table 2 of Ref. 
0, while acknowledging that most measurements are 
performed in solid matrices, which strongly influence the 
g-tensor (most notably the Ag|| components), and that 
the experimental errors are of the order of several hun- 
dred parts per million (ppm). 



Finally, by analyzing the different contributions, in- 
cluding the SOO term, we found that in all dimers apart 
from H^, the SO term accounts for more than 90% 
of 7YAg/3, and that the paramagnetic correction term 
Aggo accounts for the overwhelming majority of the SO 
term. 

To further validate our approach to the calculation of 
the g-tensor and to apply it for the first time in the solid 
state, we study two defects of a-quartz. 

The E[ center is associated with a positively charged 
oxygen vacancy, with the unpaired electron on a Si dan- 
gling bond. As in previous calculations [p|]3|,|l8|] , we model 
the defect with a 71 atom (24 Si and 47 O), positively 
charged (+1) hexagonal super-cell. We use the theoret- 
ical GGA-PBE lattice parameters (which are 1% larger 
than in experiment) and relax the atoms. For the struc- 
tural optimization we use a T only k-point sampling and 
a plane-wave cutoff of 50 Ry. The resulting relaxed struc- 
ture is very close to that of Ref. || . The EPR g-tensor is 
calculated using our relaxed structure, a plane-wave cut- 
off of 70 Ry and 4 inequivalent k-points. In Table || we 
compare our theoretical g-tensor with the experimental 
results 16 1, finding excellent agreement. 

The P2 defect center is neutral and assocated to a four- 
fold coordinated P atom, substituting for a Si atom. The 
center exists as two variants at low temperature (<140K) 
in quartz, labelled P2(I) in the ground state and P2(II) 
in the excited state pQ| . Only recently have these P de- 
fect centers been examined using DFT based total energy 
approaches |^,0|. However, up until now, the connection 
with experiment has been made using the eigenvalues of 
the hyperfine parameters alone, and the two variants of 
the defects, P2(I) and P2(II), have not been distinguished 
theoretically. 

Using the method described above, and a super-cell of 
72 atoms, 5 distinct total energy local minima are found 



as a function of the initial configuration, Table III. The 



configuration with the highest energy corresponds to a 
symmetric relaxation with the P remaining tetrahedrally 
coordinated, and O-P-0 angles of about 109°. In the 4 
other configurations the P atom moves off-center, open- 
ing up one of the 6 O-P-0 angles, which reaches a value 
of about 150°. We computed the EPR g-tensors for our 
two lowest energy structures and compare them with the 
experimental results [M in Table IV. Again, we obtain 



an excellent agreement between theory and experiment, 
which confirms that the two lowest energy theoretical 
structures correspond to the two lowest energy experi- 
mental structures. However, the comparison between g- 
tensors shows that the energy ordering between the P2(I) 
and P2(II) species is not correctly described by theory. 
This is not surprising given the small energy separation 
between the two configurations. This is expected to be 
sensitive to both the size of the super-cell and the use of 
approximated DFT functionals. 
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TABLE I. Calculated Ag tensors, in parts per million 
(ppm), for diatomic molecules. For comparison with Ref. j?J] 
we omit (in this table only) the SOO contribution to our cal- 



culations. 


A 100 Ry plane- 


-wave cutoff 


is used. 




— — ; — 

Molecule 


A S|| 




Agx 






OlrAVV 


bz FJ 


OlrAVV 


bZ, U4J 


H 2 


on 
-09 


on 
-39 


A 1 

-41 


A O 


co + 


-134 


-138 


-3223 


-3129 


CN 


-138 


-137 


-2577 


-2514 


AlO 


-141 


-142 


-2310 


-222 


BO 


-69 


-72 


-2363 


-2298 


BS 


-80 


-83 


-9901 


-9974 


MgF 


-49 


-60 


-2093 


-2178 


KrF 


-340 


-335 


61676 


60578 


XeF 


-333 


-340 


157128 


151518 



TABLE II. Calculated Ag tensors for our model E[ defect, 
and corresponding experimental data [jl6|], 



Principal values 
GIPAW (ppm) Expt. (ppm) 



Principal directions 
GIPAW Expt. 



-651 
-2255 
-2481 



-530 
-1790 
-2020 



110.0° 
142.3° 
120.4° 



223.5° 
341.6° 
121.1° 



114.5° 227.7° 
134.5° 344.4° 
125.4° 118.7° 



TABLE III. Calculated total energies, with respect to our 
lowest energy configuration. For non-tetrahedral configura- 
tions, OPO indicates the largest O-P-0 angle after the re- 
laxation, and OSiO the corresponding angle in the unrelaxed 
a-quartz structure. With 1 and s we specify whether the two 
SiO bonds forming in the OSiO angle are short or long in 
the unrelaxed structure. The number of symmetry equivalent 
structures is N c . We report our assignment of the experimen- 
tal centers based on the comparison of Table fv\. 



SiO OSiO 



OPO 



N c Energy (meV) species 



tetrahedral configuration 



I s 108.8° 
s s 108.9° 

I I 109.5° 
1 s 110.4° 



147.4° 
155.4° 
157.6° 
155.1° 



730 
243 
178 
79 




P2(I) 
P2(II) 



To summarize, we have calculated the EPR g-tensor for 
a paramagnetic defect in an extended solid for the first 
time and find our results to be in excellent agreement 
with experimental results for the E[ defect. On applying 
the method to the P2 defect, we show that comparison 
with experimental g-tensors can provide structural in- 
formation where the accuracy of DFT for energetics is 
insufficient. Combined with the calculation of hyperfine 
parameters , we expect that our GIPAW based first 
principles approach to the prediction of EPR g-tensors 
will be of great use in the assessment of models proposed 
for less well characterized paramagnetic defects, and add 
significantly to the tools available to the electronic struc- 
ture community. 

C.J. P. would like to thank the Universite Paris 6 and 
the Universite Paris 7 for support during his stay in Paris. 
The calculations were performed at the IDRIS super- 
computing center of the CNRS and on Hodgkin (SGI 
Origin 2000) at the University of Cambridge's High Per- 
formance Computing Facility. 



TABLE IV. Calculated Ag tensors for our model P2 defect 
centers, and corresponding experimental data |2o| . 

Principal values Principal directions 

GIPAW (ppm) Expt. (ppm) GIPAW Expt. 



e <t> e <j> 

Configuration P2(I) 

1249 900 65.0° 90.0° 64.8 ° 90.0° 

-980 -1100 90.0° 0.0° 90.0° 0.0° 

-3414 -3200 25.0° 270.0° 25.2° 270.0 

Configuration P2(II) 

1146 1100 47.1° 13.1° 46.5 ° 11.8° 

-824 -1000 99.0° 94.7° 101.1° 91.1° 

-3454 -3200 135.7° 355.4° 134.4° 350.1 
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